ELSEVIER 


Available  online  at  www.sciencedirect.com 


SCIENCE 


DIRECT® 


Ocean  Modelling  10  (2005)  342-368 


Ocean 

Modelling 


www.elsevier.com/locate/ ocemod 


Comparison  of  gravity  current  mixing  parameterizations 
and  calibration  using  a  high-resolution  3D  nonhydrostatic 

spectral  element  model 

Yeon  S.  Chang  a,  Xiaobiao  Xu  a,  Tamay  M.  Ozgokmen  a’*,  Eric  P.  Chassignet  a, 

Hartmut  Peters  a,  Paul  F.  Fischer  b 

a  RSMAS/MPO,  University  of  Miami,  Miami,  FL,  United  States 
b  Mathematics  and  Computer  Science  Division,  Argonne  National  Laboratory,  Argonne,  IL,  United  States 

Received  13  July  2004;  received  in  revised  form  8  November  2004;  accepted  8  November  2004 

Available  online  10  December  2004 


Abstract 

In  light  of  the  pressing  need  for  development  and  testing  of  reliable  parameterizations  of  gravity  current 
entrainment  in  ocean  general  circulation  models,  two  existing  entrainment  parameterization  schemes,  K- 
profile  parameterization  (KPP)  and  one  based  on  Turner’s  work  (TP),  are  compared  using  idealized  experi¬ 
ments  of  dense  water  flow  over  a  constant-slope  wedge  using  the  HYbrid  Coordinate  Ocean  Model 
(HYCOM).  It  is  found  that  the  gravity  current  entrainment  resulting  from  KPP  and  TP  differ  significantly 
from  one  another.  Parameters  of  KPP  and  TP  are  then  calibrated  using  results  from  the  high-order  nonhy¬ 
drostatic  spectral  element  model  Nek5000.  It  is  shown  that  a  very  good  agreement  can  be  reached  between 
the  HYCOM  simulations  with  KPP  and  TP,  even  though  these  schemes  are  quite  different  from  each  other. 
©  2004  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

Most  deep  and  intermediate  water  masses  of  the  world  ocean  are  released  into  the  large-scale 
circulation  from  high-latitude  and  marginal  seas  in  the  form  of  overflows.  For  reasons  of  mass 
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conservation,  this  downward  transport  implies  upwelling  elsewhere  in  the  ocean,  and  the  resul¬ 
ting  overturning  circulation  affects  the  large-scale  horizontal  flow  through  the  stretching  term  in 
the  vorticity  balance  (e.g.,  Gargett,  1984).  Model  representations  of  overflows  thus  determine 
more  than  just  the  properties  of  intermediate  and  deep  water  masses  in  the  ocean.  With  this 
background,  it  is  easy  to  comprehend  why  ocean  general  circulation  models  (OGCMs)  are 
highly  sensitive  to  detail  of  the  representation  of  overflows  in  these  models  (e.g.,  Willebrand 
et  al.,  2001).  Specifically,  the  entrainment  of  ambient  waters  into  overflows  is  a  prominent  oce¬ 
anic  processes  with  significant  impact  on  the  ocean  general  circulation,  and  the  climate  in 
general. 

Parameterizing  the  gravity  current  entrainment  in  coarse-resolution  OGCMs  has  proven  to  be 
challenging.  Recent  simulations  of  the  Mediterranean  overflow  employing  isopycnic  coordinates 
(Papadakis  et  al.,  2003)  and  terrain-following  coordinates  (Jungclaus  and  Mellor,  2000)  appear 
promising,  while  the  representation  of  continuous  slopes  as  steps  in  geopotential  vertical  coordi¬ 
nate  models  remains  a  daunting  problem  (e.g.,  Beckmann  and  Doscher,  1997;  Winton  et  al., 
1998;  Killworth  and  Edwards,  1999;  Nakano  and  Suginohara,  2002).  In  this  paper  we  exclusively 
focus  on  entrainment  parameterization  in  isopycnic  coordinate  models.  Isopycnic  models  have  a 
vertical  resolution  that  naturally  migrates  to  the  density  front  atop  the  gravity  current,  and  the 
amount  of  diapycnal  mixing  can  be  exactly  prescribed  (i.e.,  no  numerically-induced  diapycnal 
mixing  takes  place  as  in  geopotential  coordinate  models,  e.g.,  Griffies  et  al.,  2000).  We  conduct 
and  analyze  a  series  of  numerical  simulations  of  overflows  by  employing  parameterizations  of 
the  entrainment  simple  enough  to  be  used  in  coarse-resolution  climate  models  integrated  over 
long  time  periods.  Our  choice  of  parameterization  is  pragmatic,  motivated  by  frequent  current 
use  in  OGCMs.  Wanting  to  examine  simple  parameterizations  first,  we  deliberately  ignore  the 
more  complex,  and  computationally-expensive  schemes,  such  as  two-equation  turbulence  clo¬ 
sures,  applications  to  overflows  being  those  of  Jungclaus  and  Mellor  (2000)  and  Ezer  and  Mellor 
(2004). 

One  of  the  schemes  examined  herein  is  the  K-profile  parameterization,  KPP  (Large  et  al.,  1994, 
1999).  Its  nonlocal  treatment  of  convection  plays  no  role  in  overflows,  and  for  our  purposes,  KPP 
is  basically  a  modification  of  the  recipes  of  Pacanowski  and  Philander  (1981)  and  Munk  and 
Anderson  (1948),  the  latter  ultimately  going  back  to  observations  taken  about  a  century  ago 
and  analyzed  by  Jacobsen  (1913).  In  these  recipes,  eddy  viscosity  and  eddy  diffusivity  are  specified 
as  a  dimensional  constant  times  a  simple  analytical  function  of  the  gradient  Richardson  number. 
The  constant  is  the  maximum  possible  eddy  coefficient.  Accordingly,  the  scheme  cannot  possibly 
be  universally  valid. 

The  other  parameterization,  henceforth  referred  to  as  TP,  was  adopted  by  Hallberg  (2000)  from 
the  laboratory  experiments  of  Turner  (1986)  and  Ellison  and  Turner  (1959).  Their  original  work 
contains  an  analysis  of  the  entrainment  velocity  into  gravity  currents  as  function  of  the  bulk  Rich¬ 
ardson  number  of  this  current.  Ingeneously,  Hallberg  simply  translated  the  bulk  Richardson  num¬ 
ber  into  a  gradient  Richardson  number  ( Ri).  Rather  than  prescribing  eddy  coefficients  as  in  KPP, 
TP  thus  prescribes  the  net  entrainment  velocity  into  a  layer  as  the  velocity  difference  across  the 
layer  times  an  analytical  function  of  the  gradient  Richardson  number.  Unlike  KPP,  TP  is  hence 
proportional  to  the  forcing  by  the  shear.  TP  has  been  implemented  and  tested  in  two  isopycnic 
OGCMs,  HIM  (Hallberg  Isopycnic  coordinate  ocean  Model;  Hallberg,  2000)  and  MICOM  (Mia¬ 
mi  Isopycnic  Coordinate  Ocean  Model;  Papadakis  et  al.,  2003). 
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The  evaluation  of  the  realism  of  mixing  parameterizations  in  OGCMs  obviously  requires 
some  ground  truth.  In  this  paper,  we  bypass  the  commonly  significant  difficulties  of  comparing 
models  to  field  observations  by  taking  the  recent  three-dimensional  (3D)  high-resolution  nonhy¬ 
drostatic  simulations  of  a  generic  overflow  by  Ozgokmen  et  al.  (2004a)  as  our  ground  truth.  This 
model  resolves  the  largest  turbulent  eddies,  and,  being  nonhydrostatic,  is  physically  quite  com¬ 
plete.  Results  from  this  high-resolution,  nonhydrostatic  simulations  are  compared  to  those  from 
a  hydrostatic,  layered  OGCM  such  that  the  validity  of  the  parameterization  schemes  can  be 
examined.  Our  approach  is  as  follows:  by  comparing  the  results  from  nonhydrostatic  model 
to  those  from  OGCM,  we  quantify  the  differences  and  limitations  of  the  two  examined  Richard¬ 
son  number-dependent  parameterizations,  understand  why  and  how  these  parameterizations  can 
be  modified  to  produce  consistent  results.  Finally,  we  discuss  remaining  problems  with  both 
schemes. 

The  paper  is  organized  as  follows.  Relevant  background  information  is  given  in  Section  2.  In 
Section  3,  the  details  of  the  mixing  parameterizations  KPP  and  TP  are  introduced.  The  nonhydro¬ 
static  model  and  the  hydrostatic  OGCM  are  introduced  in  Section  4  along  with  the  experimental 
setup,  and  the  model  parameters  are  discussed.  The  results  are  presented  in  Section  5.  Finally,  the 
principal  findings  are  discussed,  and  future  directions  are  summarized  in  Section  6. 


2.  Background 

A  few  additional  remarks  on  the  physics  of  overflows  and  their  past  analyses  as  well  as  on  the 
models  employed  herein  facilitate  the  understanding  of  this  paper.  The  seminal  investigations  by 
Price  et  al.  (1993)  and  Price  and  Baringer  (1994)  reveal  that  the  mixing  of  overflows  with  the 
ambient  fluid  takes  place  over  very  small  spatial  and  time  scales.  Results  from  observational  pro¬ 
grams  in  the  Mediterranean  Sea  overflow  (Baringer  and  Price,  1997a,b),  Denmark  Strait  overflow 
(Girton  et  al.,  2001,  2003),  Red  Sea  overflow  (Peters  et  al.,  in  press;  Peters  and  Johns,  in  press), 
Faroe  Bank  Channel  (Price,  2004)  and  Antarctic  Ocean  (Gordon  et  al.,  2004)  demonstrate  the 
importance  of  small-scale  mixing  processes  in  the  dynamics  of  overflows,  and  frequently  show 
a  high  variability  of  overflow  properties  in  time  and  space.  Detailed,  quantitative  field  observa¬ 
tions  of  the  turbulent  mixing  in  overflows  are  still  few  (Johnson  et  al.,  1994a, b;  Peters  and  Johns, 
in  press). 

Hence,  much  of  our  present  understanding  of  such  mixing  is  derived  from  laboratory  tank 
experiments  (Ellison  and  Turner,  1959;  Simpson,  1969,  1982,  1987;  Britter  and  Linden,  1980; 
Turner,  1986;  Hallworth  et  al.,  1996;  Monaghan  et  al.,  1999;  Baines,  2001;  Cenedese  et  al., 
2004).  However,  when  configured  for  the  small  slopes  of  observed  overflows  [<2°],  the  dense 
source  fluid  cannot  accelerate  enough  within  the  bounds  of  typical  laboratory  tanks  [0(1  m)]  to 
produce  turbulent  behavior.  For  turbulence  to  occur,  laboratory  experiments  are  configured  with 
slopes  greater  than  10°.  It  is  further  difficult  to  maintain  a  complex  ambient  stratification  in  these 
tanks.  Ellison  and  Turner  (1959)  and  Turner  (1986)  parameterized  the  entrainment  rates  observed 
in  their  tank  experiments  as  functions  of  bulk  Richardson  numbers  of  the  flow.  Their  approach 
formed  the  basis  for  Hallberg’s  (2000)  TP  parameterization.  The  original  Turner  parameterization 
has  been  employed  in  so-called  stream  tube  models,  which  have  proven  to  be  useful  in  examining 
the  path  and  bulk  properties  of  the  Denmark  Strait  overflow  (e.g.,  Smith,  1975),  Weddell  Sea 
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overflow  (Killworth,  1977),  the  Mediterranean  overflow  (Baringer  and  Price,  1997b)  and  Red  Sea 
and  Persian  Gulf  overflows  (Bower  et  al.,  2000). 

With  the  recent  advances  in  numerical  techniques  and  computer  power,  numerical  modeling 
provides  an  alternative  avenue  to  investigate  these  processes.  Nonhydro  static,  high-resolution, 
two-dimensional  simulations  of  bottom  gravity  currents  conducted  by  Ozgokmen  and  Chassignet 
(2002)  capture  explicitly  the  major  features  of  these  currents  seen  in  laboratory  experiments,  such 
as  the  presence  of  a  head  in  the  leading  edge  and  Kelvin-Helmholtz  vortices  in  the  trailing  fluid. 
Subsequently,  this  model  was  used  to  simulate  the  part  of  the  Red  Sea  outflow  in  a  submarine 
canyon,  which  naturally  restricts  motion  in  the  lateral  direction  such  that  the  use  of  a  two-dimen¬ 
sional  (2D)  model  provides  a  reasonable  approximation  to  the  dynamics.  It  was  shown  (Ozgok¬ 
men  et  al.,  2003)  that  this  model  adequately  captures  the  general  characteristics  of  mixing  in  the 
Red  Sea  overflow  within  the  limitations  of  a  2D  model.  These  limitations  include  lack  of  edge  ef¬ 
fects  or  intrusions  from  channel  walls  associated  with  the  spanwise  structure.  Recently,  a  parallel 
high-order  spectral  element  Navier-Stokes  solver,  Nek5000,  developed  by  Fischer  (1997),  was 
used  to  conduct  nonhydrostatic  3D  simulations  of  bottom  gravity  currents  (Ozgokmen  et  al., 


2004a, b). 


In  this  study,  our  objective  is  to  explore  how  mixing  parameterizations  perform  in  an  idealized 
setting  that  represents  the  very  basics  of  shear-induced  mixing  in  bottom  gravity  currents,  e.g. 
flow  of  a  dense  water  mass  released  at  the  top  of  a  sloping  wedge.  To  this  end,  we  conduct  exper¬ 
iments  with  a  layered  hydrostatic  OGCM,  HYCOM  (HYbrid  Coordinate  Ocean  Model),  employ¬ 
ing  KPP  and  TP,  and  compare  the  results  to  those  from  high-resolution  3D  nonhydrostatic 
simulations  by  Ozgokmen  et  al.  (2004a).  We  start  out  from  the  questions  (a)  how  the  results  from 
KPP  and  TP  compare  to  each  other  and  to  those  from  the  3D  nonhydrostatic  simulations,  and  (b) 
how  the  results  change  and/or  converge  as  a  function  of  the  grid  spacing.  Discussing  our  results, 
we  examine  what  the  principal  limitations  of  the  KPP  and  TP  parameterizations  are  and  how  can 
they  be  developed  into  a  consistent  formulation  for  use  in  layered  ocean  models. 


3.  Mixing  parameterizations 


3.1.  KPP 


The  K-profile  parameterization  (Large  et  al.,  1994,  1999;  KPP)  provides  a  prescription  for  mix¬ 
ing  from  surface  to  bottom,  smoothly  matching  large  values  of  the  eddy  diffusivity  ( K)  in  the  sur¬ 
face  boundary  layer  to  small  values  in  the  interior  of  the  ocean.  KPP  has  been  popular  because  it 
includes  prescriptions  for  a  fairly  wide  range  of  physical  processes,  shear-driven  mixing  in  low -Ri 
regions,  constant  background  internal- waves  induced  mixing  allowing  counter-gradient  fluxes. 
Herein,  only  the  shear-induced  mixing  is  important,  a  component  of  KPP  that  was  not  specifically 
tailored  to  gravity  currents.  Shear-driven  mixing  is  expressed  in  terms  of  the  gradient  Richardson 
number  Ri  calculated  at  layer  interfaces: 


Ri  =  N2 


(1) 
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where  the  numerator  is  the  buoyancy  frequency,  N2  =  —  ^  written  here  for  an  incompressible 
fluid  for  simplicity,  and  where  g  =  9.81  m2  s  1  is  the  gravitational  acceleration.  The  denominator 
in  (1)  is  the  vertical  shear.  The  vertical  diffusivity  is  then  related  to  Ri  as 

Ri 
Ric 

so  that  vertical  diffusivity  is  zero  when  Ri  5s  Ric  corresponding  to  the  case  in  which  stratification 
overcomes  the  effect  of  vertical  shear  and  prohibits  vertical  mixing  (Fig.  la).  In  this  case,  mixing 
can  only  take  place  in  the  horizontal  plane,  or  as  so-called  “pancake”  mixing  demonstrated  in 
laboratory  experiments  of  stratified  flow  (e.g.,  Fernando,  2000).  Ric  is  set  to  Ric  =  0.7. 

With  decreasing  Ri  <  Ric,  the  vertical  diffusivity  coefficient  gradually  increases  to  account  for 
mixing  induced  by  high  vertical  shear  and/or  weaker  stratification.  At  the  limit  of  Ri  =  0,  mixing 
takes  place  as  in  homogeneous  (unstratified)  fluid.  The  concept  behind  this  component  of  KPP  is 
taken  from  Pacanowski  and  Philander  (1981),  but  the  shape  of  the  mixing  curve  in  KPP  was  ad¬ 
justed  to  show  better  agreement  with  observational  results  from  equatorial  mixing  by  Gregg  et  al. 


K, 


shear 


=  KU 


1  —  min  (  1 , 


KPP 


TP 


Fig.  1 .  Mixing  curves  of  KPP  and  TP.  (a)  ^Shear  (cm2  s  *)  as  a  function  of  Ri  in  KPP,  and  (b)  we/A  U  as  a  function  of  Ri 
in  TP. 
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(1985)  and  Peters  and  Gregg  (1988)  for  the  regime  of  0.3  ^  Ri  ^  0.7.  The  diffusivity  value  at 
Ri  =  0  was  determined  from  large  eddy  simulation  (LES)  studies  of  the  upper  tropical  ocean  as 
Kmax  =  50  cm2  s-1  (e.g.,  Large,  1998). 

The  boundary  shear  stress  component  of  KPP  can  be  invoked  to  account  for  the  bottom  stress 
in  the  bottom  boundary  layer  (BBL): 

■^shear  =  ^b^sGs.  (3) 

Here,  ///„  is  the  bottom  boundary  layer  thickness,  estimated  as  the  total  thickness  of  layers, 
counted  from  the  bottom,  in  which  the  Richardson  number  remains  lower  than  the  critical  value. 
The  velocity  scale  u  s  is  a  linear  function  of  the  friction  velocity  which  is  proportional  to  the 
square  root  of  the  bottom  stress  pcd|ut>|ub,  where  =  3  x  10-3  and  Ub  is  the  bottom  velocity. 
G$  is  a  third-order  smooth  shape  function  to  match  with  the  interior  profiles.  The  reader  is  re¬ 
ferred  to  Large  et  al.  (1994)  for  further  detail  on  KPP,  and  to  Halliwell  (2004)  on  the  numerical 
implementation  of  KPP  in  HYCOM. 

3.2.  TP 

Based  on  the  results  and  parameterization  by  Turner  (1986)  and  Ellison  and  Turner  (1959)  and 
Hallberg  (2000)  developed  a  mixing  parameterization,  in  which  the  net  entrainment  velocity  into 
layers  in  gravity  currents  wE  is  expressed  as 

f  CA  At/(0.08  —  0.1 /?!')/ (1  +  5 i?/)  if  Ri  <  0.8, 

WE=\o  if  Ri  >0.8,  (4) 

where  AU  is  the  mean  velocity  difference  across  layers,  CA  =  1  is  a  proportionality  constant,  and 
the  cut-off  Richardson  number  is  Ric  =  0.8.  The  reader  is  referred  to  Hallberg  (2000)  for  further 
detail  of  the  numerical  implementation. 

Because  both  TP  and  KPP  employ  functions  of  Ri,  there  is  similarity  between  the  prescriptions 
of  shear-driven  mixing  in  TP  (4)  and  KPP  (2).  Fig.  la  and  b  depict  the  different  shapes  of  the  mix¬ 
ing  curves  and  the  different  values  of  the  critical  Ri.  An  important  difference  lies  in  the  hard  limit 
K  ^  KmSiX  m  KPP  and  the  absence  of  any  such  limit  m  TP . 


4.  The  numerical  models  and  experimental  configuration 

4.1.  Nonhydrostatic  3D  model  Nek5000 

Results  from  the  high-order  parallel  spectral  element  Navier-Stokes  solver,  Nek5000,  are  used 
as  a  reference.  This  model  is  documented  in  detail  by  Fischer  (1997),  Fischer  et  al.  (2000),  Tufo 
and  Fischer  (1999)  and  Fischer  and  Mullen  (2001).  A  short  description  of  the  model  in  the  context 
of  bottom  gravity  current  experiments  can  also  be  found  in  Ozgokmen  et  al.  (2004a). 

Nek5000  is  a  state-of-the-art  general  computational  fluid  dynamics  code  (see  http://www-unix. 
mcs.anl.gov/~fischer/  for  applications)  that  is  based  on  the  spectral  element  method  (SEM),  and 
offers  several  fundamental  advantages  with  respect  to  the  more  common  numerical  discretization 
techniques  (finite-difference,  finite-element,  finite-volume  and  spectral);  (a)  SEM  combines  the 
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geometrical  flexibility  of  finite  element  method  with  the  numerical  accuracy  of  spectral  expansion 
(e.g.,  the  geometrical  flexibility  of  SEM  has  been  exploited  to  explore  the  behavior  of  bottom 
gravity  currents  over  complex  topography  in  Ozgokmen  et  al.,  2004b);  (b)  SEM  has  minimal  dis¬ 
sipation  and  dispersion  errors,  which  are  important  in  problems  involving  propagation  of  high 
gradients  and  mixing,  such  as  in  the  present  problem;  (c)  SEM  provides  dual  path  to  convergence, 
either  via  elemental  grid  refinement  or  by  increasing  the  polynomial  degree;  (d)  SEM  offers  com¬ 
putational  advantages  for  scalability  on  parallel  computers  (Tufo  and  Fischer,  1999). 

In  the  present  setup,  Nek5000  is  configured  to  solve  the  Boussinesq  equations: 

^ =-VP  +  Vlu-RaSz ,  (5) 

V  •  u  =  0,  (6) 


D  S 

d7 


Pr~x  V2rS, 


where  the  material  (total)  derivative  is  ^ 

2  S'  S'  92 

V2  := - 1 - h  r — . 

r  dx2  dy2  0z2 


^  +  u  •  V,  and  the  anisotropic  diffusivity  is 


(7) 

(8) 


The  variables  are  the  velocity  vector  field  u  =  («,  v,  w)  and  the  pressure  p,  and  the  nondimensional 
parameters  are  Ra  =  ( gfiASH3)/v2h  the  Rayleigh  number,  the  ratio  of  the  strengths  of  buoyancy 
and  viscous  forces,  where  H  is  the  domain  depth  and  AS  is  the  salinity  range  in  the  system, 
P  =  7  x  10-4  psu-1  is  the  salinity  contraction  coefficient  (a  linear  equation  of  state  is  used); 
Pr  =  vh/Kh  the  Prandtl  number,  the  ratio  of  viscous  and  saline  diffusion;  and  r  =  vv/vh  =  KJK h, 
the  ratio  of  vertical  and  horizontal  diffusivities. 


4.2.  Hydrostatic  ocean  model  HYCOM 

The  development  of  Hybrid  Coordinate  Ocean  Model  (HYCOM)  was  motivated  by  the  fact 
that  no  single  vertical  coordinate — depth,  density,  or  terrain-following — can  be  by  itself  optimal 
everywhere  in  the  ocean.  The  default  configuration  in  HYCOM  is  one  that  is  isopycnal  in  the 
open  stratified  ocean  but  smoothly  reverts  to  a  terrain-following  coordinate  in  shallow  coastal  re¬ 
gions  and  to  fixed  pressure-level  coordinates  in  the  surface  mixed  layer  and  unstratified  seas.  In 
doing  so,  the  model  ideally  combines  the  advantages  of  the  different  types  of  coordinates  in  sim¬ 
ulating  coastal  and  open  ocean  circulation  features.  The  basic  principles  of  this  generalized  ver¬ 
tical  coordinate  model  are  described  in  Bleck  (2002),  Chassignet  et  al.  (2003)  and  Halliwell 
(2004),  and  detailed  documentation  is  readily  available  from  http://hycom.rsmas.miami.edu. 

4.3.  Experimental  configuration 

In  Nek5000,  the  model  domain  has  a  horizontal  (streamwise)  length  of  Lx-  10  km  and  a  span- 
wise  width  of  Ly  =  2  km.  The  depth  of  the  water  column  ranges  from  400  m  at  x  =  0  to 
H  =  1000  m  at  x  =  10  km,  hence  the  background  slope  angle  is  9  =  3.5°  (Fig.  2a).  The  boundary 
conditions  at  the  bottom  are  no-slip  and  no-normal  flow  for  the  velocity  components,  and  no- 


Y.S.  Chang  et  al.  /  Ocean  Modelling  10  (2005)  342-368 


349 


Dirichlet 


*'/0 
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Fig.  2.  Configuration  of  experiments  in  Nek5000.  (a)  Schematic  depiction  of  the  domain  geometry  and  boundary 
conditions  (length  scale  is  in  km),  and  (b)  velocity  profile  at  the  forcing  boundary  and  the  initial  distribution  of  salinity. 
Distribution  of  elements  is  depicted  in  the  background. 


normal  flux,  0570 n  =  0,  with  n  being  the  normal  direction  to  the  boundary,  for  salinity.  Rigid-lid 
and  free-slip  boundary  conditions  are  used  at  the  top  boundary.  The  model  is  driven  by  the  veloc¬ 
ity  and  salinity  forcing  profiles  at  the  inlet  boundary  at  v  =  0.  The  model  is  initialized  by  using  a 
salinity  distribution  of  the  (dimensional)  form 

o>=f  «p(-{M1+0;rtn(,a)}  J 


1  —  cos 


H 


0.4// 
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The  sinusoidal  perturbation  in  the  spanwise  direction  facilitates  the  transition  to  3D  flow.  A  sec¬ 
tion  of  the  initial  salinity  profile  across  the  middle  of  the.  domain  (Fig.  2b)  shows  that  the  initial 
thickness  of  the  dense  water  mass  is  h0  =  200  m.  The  velocity  distribution  at  the  inlet  matches  no¬ 
slip  at  the  bottom  and  free-slip  at  the  top  in  such  a  way  that  the  depth-integrated  mass  flux  across 
this  boundary  is  zero.  The  amplitude  of  the  inflow  velocity  profile  scales  with  the  propagation 
speed  of  the  gravity  current,  which  is  zero  at  t  -  0  and  reaches  a  constant  value  UF  shortly  after 
its  release  such  that  the  bulk  Froude  number  is  near  critical,  Fr  =  UF/^/g[3ASho  «  0.9,  which  is 
characteristic  for  overflows  emanating  from  narrow  straits  (e.g.,  Price  and  Baringer,  1994;  Murray 
and  Johns,  1997).  As  the  interior  is  initialized  with  homogeneous,  light  ( S  =  0)  water,  the  density 
front  propagation  is  the  fastest  signal  in  the  system.  Density  currents  reach  the  exit  boundary  after 
about  10,000  s,  at  which  point  the  integrations  are  terminated  such  that  potential  complications 
involving  the  outflow  boundary  are  avoided,  albeit  with  the  limitation  of  focusing  only  on  the 
start-up  phase  of  plumes  rather  than  those  in  a  statistically  steady  state.  Finally,  periodic  bound¬ 
ary  conditions  are  applied  at  the  lateral  boundary.  The  domain  is  discretized  using  4000  elements 
with  lOth-order  polynomials  in  each  spatial  direction  within  the  elements,  hence  a  total  of  4  x  106 
grid  points  are  employed.  The  remaining  model  parameters  are  listed  in  Table  1  and  the  reader  is 
referred  to  Ozgokmen  et  al.  (2004a)  for  more  detailed  discussion.  The  calculations  were  carried 
out  on  a  Linux  cluster  running  on  32  Athlon  1.7  GHz  processors,  and  it  takes  approximately  9 
days  to  complete  (simulated  to  real-time  ratio  of  ^1/60). 

In  HYCOM  experiments,  the  model  parameters  and  the  physical  conditions  of  the  model  are 
set  as  closely  as  possible  to  those  in  Nek5000  for  the  comparison.  Two  major  changes  are  made 
in  the  domain  configuration  of  HYCOM  experiments.  First,  the  sloping  portion  of  the  domain  is 
extended  from  10  km  to  20  km  (while  keeping  6  =  3.5°)  in  order  to  obtain  more  reliable  estimates 
of  the  entrainment  parameters  (Fig.  3a).  Second,  an  inlet  system  is  designed  using  a  reservoir  of 
10  km  length,  initially  filled  with  dense  water  (Fig.  3b).  Shortly  after  the  dam  break  (at  t  -  15  min, 
Fig.  3c),  the  system  becomes  nearly  equivalent  to  the  initial  conditions  used  in  Nek5000  (Fig.  lb) 
without  the  use  of  velocity  boundary  conditions,  and  in  a  fashion  consistent  with  changes  in  grid 
spacing.  The  reservoir  contains  an  adequate  amount  of  dense  water  for  the  sloping  portion  of  the 
domain.  Free  surface  boundary  conditions  are  used  at  the  top  boundary,  and  a  quadratic  drag  law 
with  a  drag  coefficient  of  cd  =  3  x  10-3  is  applied  at  the  bottom.  Free-slip  boundary  conditions  are 
applied  at  the  lateral  boundaries. 


Table  1 


Parameters  of  the  Nek5000  nonhydro  static  model  simulation 


Domain  Size  (Lx,  Lz  -  H ,  Ly) 

(104  m,  103  m,  2  x  103  m) 

Slope  angle 

0  =  3.5° 

Rayleigh  number 

Ra  =  5  x  106 

Prandtl  number 

Pr=  1 

Diffusivity  ratio 

r  =  2  x  10-2 

Salinity  range 

AS  =  1.0  psu 

Number  of  elements  (x,z,y) 

50,  8,  10 

Polynomial  degree 

N  =  10 

Number  of  grid  points 

4x  106 

Time  step 

At  =  0.85  s 

Y.S.  Chang  et  al.  /  Ocean  Modelling  10  (2005)  342-368 


351 


N  -| 

(b)  1.5 


x,  km 


Fig.  3.  Computational  domain  and  initial  conditions  in  HYCOM  experiments,  (a)  Domain  geometry  in  3D,  (b)  the 
initial  salinity  distribution  (x-z  plane),  and  (c)  salinity  distribution  t  =  15  min. 


One  of  the  objectives  of  the  present  study  is  to  explore  the  behavior  of  gravity  current  simula¬ 
tions  in  HYCOM  as  a  function  of  grid  spacing.  Five  different  horizontal  resolutions  are  used: 
Ax  =  Ay  =  1000,500, 100,50,  and  20  m.  The  gradual  strengthening  of  nonhydrostatic  effects  limi¬ 
ted  the  minimum  grid  spacing  to  Ax  =  20  m.  The  horizontal  viscosity  changes  in  HYCOM  in  pro¬ 
portion  to  grid  spacing  as  v  =  max{wd  Ax,  [(^  —  |})2  +  +  ^)2]1^2  Ax2},  where  u^-2  cm  s_1. 

Finally,  experiments  were  conducted  for  both  5  ana  1 1  layers,  but  since  there  was  no  significant 
difference,  we  only  present  results  from  the  5-layer  runs  in  the  following.  The  main  parameters  for 
HYCOM  experiments  are  summarized  in  Table  2. 

One  important  factor  in  the  dynamics  of  oceanic  overflows  is  rotation.  The  scale  at  which  the 
Coriolis  force  becomes  comparable  to  the  buoyancy  force  is  a  complex  function  of  the  slope  angle, 
stratification,  and  friction  (e.g.,  Griffiths,  1986).  A  simple  spatial  scale  for  rotational  effects  to  be¬ 
come  important  is  given  by  the  radius  of  deformation  y/gH/ f,  which,  with  the  stated  experimental 
parameters,  is  approximately  17  km  at  midlatitude,  as  compared  to  domain  lengths  of  10  or 
20  km.  The  rotation  time  scale  is  f^1  ~  15,000  s,  while,  as  shown  below,  the  gravity  currents  cross 
the  domain  in  ^10,000  or  20,000  s.  Hence,  the  results  presented  here  apply  to  the  phase  before  the 
impact  of  rotation  starts  influencing  the  flow  patterns. 
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Table  2 

Parameters  of  the  HYCOM  simulations 


KPP 

TP 

Domain  size  ( Lx ,  Lz  =  H ,  Ly) 

(3  x  104  m,  1.6  x  103  m,  2  x  103  m) 

Same 

Slope  angle 

0  =  3.5° 

Same 

Mixing  parameters 

Imax  -  50,  2500  cm2  s-1 

Ca  =  1,  0.15 

Salinity  range 

AS  =  1.0  psu 

Same 

Horizontal  resolution 

Ax  =  Ay  =  20,  50,  100,  500,  1000  m 

Same 

Vertical  resolution 

5,  11  layers 

Same 

Time  step 

At=  1,  9  s 

Same 

Given  that  in  some  overflows  the  bulk  of  the  entrainment  takes  place  over  a  very  small  distance 
(e.g.,  over  a  distance  of  approximately  40-50  km  in  the  Mediterranean  Sea  over-flow  according  to 
Fig.  7a  of  Baringer  and  Price,  1997b),  in  the  high-resolution  studies  of  [49-52]  it  was  considered 
important  that  the  detail  of  such  entrainment  be  captured.  In  this  sense  this  study  complements 
other  process  studies  focusing  on  the  larger-scale  behavior  (e.g.,  Ezer  and  Mellor,  2004). 


5.  Results 

5.7.  Description 

The  evolution  of  the  salinity  distribution  in  the  Nek5000-experiments  is  shown  in  Fig.  4.  The 
system  is  initialized  as  described  in  Section  4.3.  The  initial  development  of  the  system  is  that  of 
the  so-called  lock-exchange  flow  (e.g.,  Keulegan,  1958;  Simpson,  1987),  in  which  the  lighter  fluid 
remains  on  top  and  the  denser  overflow  propagates  downslope.  The  dense  gravity  current  quickly 
develops  a  characteristic  “head”  at  the  leading  edge  of  the  current  (Fig.  4b).  The  head  is  half  of  a 
dipolar  vortex,  which  is  a  generic  flow  pattern  that  tends  to  form  in  two-dimensional  systems  by 
self-organization  of  the  flow  (e.g.,  Flierl  et  al.,  1981;  Nielsen  and  Rasmussen,  1996),  and  which 
corresponds  to  the  most  probable  equilibrium  state  maximizing  entropy  (Smith,  1991).  The  head 
grows  and  is  diluted  as  the  gravity  current  travels  further  down  the  slope,  the  result  of  entrainment 
of  fresh  ambient  fluid.  The  flow  along  the  leading  edge  of  the  current  is  composed  of  a  complex 
pattern  of  so-called  lobes  and  clefts  that  are  highly  unsteady  (Fig.  4c  and  d)  and  well-known  fea¬ 
tures  from  laboratory  experiments  tracing  back  to  the  work  of  Simpson  (1972).  It  was  conjec¬ 
tured,  e.g.,  by  Simpson  (1987)  that  a  gravitational  rise  of  the  thin  layer  of  light  fluid  which  the 
gravity  current  overruns  is  responsible  for  the  breakdown  of  the  flow  at  the  leading  edge.  Re¬ 
cently,  Hartel  et  al.  (2000)  put  forth  that  instability  associated  with  the  unstable  stratification  pre¬ 
vailing  at  the  leading  edge  between  the  nose  and  stagnation  point  of  the  front  could  also  account 
for  this  behavior. 

In  the  trailing  fluid,  the  initial  instabilities  appear  to  be  2D  Kelvin-Helmholtz  rolls  that  span 
the  entire  width  of  the  domain  (Fig.  4c).  These  rolls  gradually  exhibit  transition  to  3D  (Fig. 
4d).  The  development  of  spanwise  instabilities  in  Kelvin-Helmholtz  rolls  was  investigated  by  Kla- 
assen  and  Peltier  (1991),  who  classified  them  into  two  categories.  The  first  are  dynamical  second¬ 
ary  instabilities  that  tend  to  initiate  in  the  vortex  core  and  the  interface  between  strongly 


Y.S.  Chang  et  al.  /  Ocean  Modelling  10  (2005)  342-368 


353 


(a)  (b) 


Fig.  4.  Distribution  of  salinity  surface  0.3  ^  S  <  0.6  in  Nek5000  at  (a)  t  =  0,  (b)  t  -  2125  s  (^0.6  h),  (c)  t  =  4675  s 
(«1.3h),  (d)  t  =  9350  s  («2.6  h). 

rotational  and  weakly  rotational  fluid  and  that  develop  independently  at  different  growth  rates. 
The  second  category  are  convective  secondary  instabilities  in  the  statically  unstable  regions,  which 
develop  as  the  interface  between  the  two  streams  overturns.  Because  of  these  instabilities,  Kelvin- 
Helmholtz  billows  cannot  preserve  their  coherence  in  the  lateral  y-direction,  and  break  down.  In 
contrast,  Kelvin-Helmholtz  billows  in  2D  can  grow  by  pairing  (e.g.,  Corcos  and  Sherman,  1984). 
Therefore,  Kelvin-Helmholtz  billows  in  3D  result  in  smaller  coherent  overturning  structures  than 
those  in  2D,  and  consequently,  the  entrainment  parameter  in  3D  simulations  was  found  to  be 
smaller  than  that  in  2D  (Ozgokmen  et  al.,  2004a).  Finally,  a  span  wise  averaged  salinity  distribu¬ 
tion  is  depicted  in  Fig.  5  for  a  visual  comparison  to  results  from  HYCOM  experiments. 


Fig.  5.  Distribution  of  spanwise  averaged  salinity  in  Nek5000  at  t  -  9350  s. 
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In  HYCOM  experiments,  mixing  of  the  bottom  gravity  current  with  the  ambient  fluid  entirely 
depends  on  the  mixing  parameterizations,  in  the  absence  of  which,  none  would  occur  by  design. 
Results  using  KPP  and  TP  at  coarse  resolution  of  Ax  =  1000  m  are  illustrated  in  Fig.  6.  At 
t  -  4950  s  (Fig.  6a,  corresponding  approximately  to  Fig.  4c;  note  that  in  HYCOM  the  sloping  part 
of  the  domain  is  twice  as  long  as  that  in  Nek5000),  there  is  formation  of  a  characteristic  head 
using  KPP.  The  head  grows  in  time  and  breaks  into  two  parts  before  the  gravity  current  reaches 
the  end  of  the  domain  at  /  =  13,050  s  (Fig.  6b  and  c).  In  contrast,  the  growth  of  the  head  is  nearly 
absent  in  the  experiment  with  TP  (Fig.  6d-f).  Clearly,  there  is  significantly  more  entrainment  and 
dilution  of  the  gravity  current  resulting  in  a  much  slower  propagation  speed  of  the  gravity  current 
in  the  experiment  with  TP  than  that  with  KPP.  There  is  no  significant  variation  in  y-direction  in 
either  experiment  so  that  they  are  effectively  2D,  as  are  the  geometry,  forcing  and  initial 
conditions. 

Results  from  KPP  and  TP  at  fine  resolution  of  Ax  =  20  m  are  shown  in  Fig.  7.  At  this  resolu¬ 
tion,  KPP  results  in  a  great  deal  of  fine-scale  structure  in  the  tail  of  the  gravity  current.  As  shown 
by  Gallacher  and  Piacsek  (in  press),  this  fine-scale  behavior  occurs  because  of  the  hydrostatic 
approximation,  which  is  shown  to  lead  to  unphysical  noise  related  to  the  overestimation  of  the 
vertical  velocity  at  high  horizontal  resolution.  (This  noise  is  more  significant  in  KPP  than  in 
TP,  but  both  cases  are  unstable  for  Ax  ^  10  m.)  TP  at  this  resolution  yields  a  head  at  the  leading 
edge,  which  is  smaller  than  that  in  KPP.  The  main  result  remains  the  same  in  general;  there  ap¬ 
pears  to  be  a  significant  difference  in  entrainment  resulting  from  these  two  schemes  in  that  TP  re¬ 
sults  into  substantially  more  entrainment  than  KPP. 


KPP: 


TP: 


(a)  salinity  in  2-D  (X-Z),  t=4950  sec 


(d)  salinity  in  2-D  (X-Z),  t=4950  sec 


10  15  20  25 


x,  km 


x,  km 


Fig.  6.  Evolution  of  the  distribution  of  salinity  perturbation  (S'  =  S  —  So)  in  time  in  HYCOM  experiments  (a) 
t  =  4960  s,  (b)  t  =  9450  s,  and  (c)  t  =  13,050  s  using  KPP  with  Ax  =  1000  m,  and  (d),  (e),  (f)  using  TP. 


Y.S.  Chang  et  al.  /  Ocean  Modelling  10  (2005)  342-368 


355 


KPP:  TP: 
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Fig.  7.  Evolution  of  salinity  distribution  in  time  in  HYCOM  experiments  (a)  spanwise-average  at  t  -  4960  s,  (b) 
spanwise-average  at  t  =  9450  s,  and  (c)  spanwise-average  at  t  =  13,050  s  using  KPP  with  Ax  =  20  m,  and  (d),  (e),  (f) 
using  TP. 

It  is  also  well-known  that  bottom  gravity  currents  propagate  with  a  constant  speed  provided 
that  the  input  flux  is  constant.  In  this  flow  the  gravitational  force  is  balanced  by  a  combination 
of  bottom  friction  and  entrainment  (e.g.,  Britter  and  Linden,  1980).  Fig.  8  shows  the  position 
of  the  density  front  as  a  function  of  time,  XF(t),  in  experiments  with  KPP  and  TP.  In  all  exper¬ 
iments,  the  speed  of  propagation,  UF  =  dXF(t)ldt,  is  approximately  constant.  When  scaled  with 
the  speed  of  long  internal  waves,  y/g'ho  =1.18  m  s-1,  where  g'  =  gfiAS  «  7  x  10-3  m  s-2  and 
h0  =  200  m,  UF/^/g'h0  =  0.92  in  the  case-of-TP,  and  UF/^/g'h0  =  1.23  in  the  case  of  KPP  in 
coarse  resolution  (Ax  =  1000  m)  experiments.  This  difference  in  propagation  speed  further  empha¬ 
sizes  the  difference  in  total  entrainment  from  TP  and  KPP  schemes. 


5.2.  Entrainment 


Turner  (1986)  defined  an  entrainment  parameter  in  bottom  gravity  currents  as  the  change  of  the 
dense  flow  thickness  h  along  the  streamwise  direction  X, 


dh 

=  dX  ’ 

a  2D  expression,  which  can  be  mapped  to  3D  flows  as  (Ozgokmen  et  al.,  2004a), 
h(t)  -  h0(t ) 


(9) 


m  = 


m 


(10) 
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Fig.  8.  Position  of  the  salinity  front,  XF  (in  m),  as  a  function  of  time  in  HYCOM  experiments  with  (a)  Ax  =  1000  m. 
Blue  line:  KPP  with  Kmax  =  50  cm2  s_1,  green:  KPP  with  Imax  =  2500  cm2  s_1,  red:  TP,  with  CA  =  1.0,  and  black:  TP 
with  CA  =  0.15.  (For  interpretation  of  the  references  in  colour  in  this  figure  legend,  the  reader  is  referred  to  the  web 
version  of  this  article.) 


where  £(t)  =  L~l  YF(y\  t) cl/  —  x0  is  the  spanwise-averaged  length  of  the  dense  overflow  mea¬ 
sured  between  the  reference  station  x0  and  the  leading  edge  XF(y,  t),  h(t)  is  the  total  (with  entrain¬ 
ment)  mean  thickness  estimated  from 

1  r-Ly  r-XF(y',t) 

J  h{x',y',t)dx'dy',  (11) 

between  a  reference  station  of  x0  and  the  leading  edge  of  the  density  current  XF.  The  overflow 
thickness  h  is  calculated  from 


r 

h(x,y,t)=  /  d(x,y,zr  ,t)dzf  where  S(x,y,z,  t) 

Jo 


0,  when  S(x,y,z,t)  <  e, 
1,  when  S(x,y,z,t)  ^  e. 


(12) 


The  salinity  interface  threshold  value  is  taken  as  e  =  0.2  (psu),  since  in  the  case  of  Nek5000,  it  is 
the  lowest  salinity  value  remaining  as  a  coherent  part  of  the  gravity  current  (fluid  particles  with 
lower  salinity  tend  to  detach  from  the  current  and  be  advected  with  the  overlying  counter  flow). 
Finally,  h0(t)  is  the  original  (without  any  entrainment)  mean  thickness  estimated  from 

hit)  [  [  [  u(x0,y',z?,t')dz'dy'dt'.  (13) 

tytj-Ly  J o  Jo  Jz°+h 


As  E(t)  accounts  for  the  mixing  process  of  the  gravity  current  with  the  ambient  flows,  the  com¬ 
parison  of  E(t)  between  the  results  from  HYCOM  and  Nek5000  can  diagnose  the  appropriateness 
of  the  vertical  mixing  parameters  in  the  hydrostatic  ocean  models.  In  Ozgokmen  et  al.  (2004a),  it 
was  found  that  the  entrainment  parameter  converges  to  E  «  4  x  10-3  in  3D  experiments,  whereas 
E  k,  9  x  10-3  in  2D,  because  of  the  differences  between  2D  and  3D  turbulence  discussed  above. 
These  results  are  plotted  in  Fig.  9  in  comparison  to  HY COM  experiments  with  KPP  and  TP  carried 
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E(t)  in  HYCOM-KPP  vs  Nek5000 


E(t)  in  HYCOM-TP  vs  Nek5000 


Fig.  9.  Time  evolution  of  entrainment  parameters,  E(t),  in  HYCOM  experiments  with  (a)  KPP  and  (b)  TP  at  different 
horizontal  grid  spacings;  (x)  Ax  =  1000  m,  (A)  Ax  =  500  m,  (O)  Ax  =  100  m,  (0)  Ax  =  50  m,  (□)  Ax  =  20  m. 
Entrainment  parameters  from  2D  and  3D  nonhydro  static  experiments  with  Nek5000  are  shown  in  the  background, 
dotted:  2D,  dashed:  3D. 


out  using  five  different  horizontal  grid  resolutions.  Fig.  9a  illustrates  that  the  entrainment  para¬ 
meter  of  KPP  converges  to  a  mean  value  of  E  «  1  x  10-3,  which  is  a  significantly  less  than  the  value 
obtained  from  Nek5000.  This  result  is  reasonably  independent  of  the  horizontal  grid  scale;  while 
some  oscillatory  behavior  is  evident  with  Av  =  1000  m,  and  some  underestimate  is  obtained  with 
Av  =  500  m,  the  entrainment  profiles  are  very  similar  with  the  resolutions  of  100,  50,  and  20  m. 
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In  contrast,  E(t)  obtained  from  HYCOM  with  TP  converges  to  E  «  8.5  x  10-3  (Fig.  9b),  which 
is  in  very  good  agreement  with  the  estimate  given  by  Turner  as  E  -  (5  +  9)  x  10-3  for  9  =  3.5°. 
While  this  entrainment  parameter  is  in  approximate  agreement  with  the  2D  results  from 
Nek5000,  it  is  larger  by  a  factor  of  2  than  the  value  from  the  3D  simulations.  Similar  behavior 
to  KPP  with  regard  to  the  grid  scale  applies  to  TP  as  well;  oscillatory  convergence  is  displayed 
with  Ax  =  1000  m,  slight  underestimation  of  E  results  from  Ax  =  500  m,  and  nearly-identical  re¬ 
sults  are  obtained  with  Ax  =  100,  50,  and  20  m.  By  conducting  experiments  both  with  At  =  9  s  and 
with  At  =  1  s,  it  is  found  that  the  oscillatory  convergence  of  E(t)  at  Ax  =  1000  m  is  entirely  related 
to  spatial  discretization. 

5.3.  Calibration  of  KPP  and  TP  using  Nek5000 

The  computations  discussed  in  Section  5.2  clearly  illustrate  that,  with  their  original  parameters, 
KPP  and  TP  lead  to  significantly  different  results.  Further  evidence  why  these  parameters  should 
be  considered  tunable  and  why  there  is  merit  in  considering  a  calibration  of  the  parameters  of 
both  KPP  and  TP  relative  to  results  from  Nek5000  runs  is  provided  in  the  summary  section  fur¬ 
ther  below. 

Both  KPP  and  TP  have  two  parameters,  the  critical  gradient  Richardson  number,  above  which 
turbulent  mixing  terminates,  and  an  amplitude  parameter.  We  have  not  found  significant  devia¬ 
tions  from  the  results  shown  in  Fig.  9a  and  b  when  Ric  is  changed  by  ±50%.  However,  varying  the 
respective  amplitude  parameters  in  TP  and  KPP  causes  large  variations  in  model  outputs.  Fig.  10 
shows  the  results  obtained  by  adjusting  the  coefficients  Km.dX  and  Ca  in  KPP  and  TP,  respectively, 
such  that  E(t)  approximately  matches  that  obtained  from  3D  Nek5000  experiments.  Specifically, 
the  rms  deviation  between  E{t)  in  HYCOM  with  resolutions  of  Ax  =  100  m,  50  m,  20  m  is  mini¬ 
mized  for  7000  s  <  /  <  10,000  s.  This  calibration  results  in  the  optimized  coefficient  values  of 
Xmax  -  2500  cm2  s-1  for  KPP  and  Ca  =  0.15  for  TP.  Using  these  parameters,  snapshots  of  the 
gravity  current  at  the  different  times,  and  from  different  horizontal  resolutions,  show  good  agree¬ 
ment  between  the  results  obtained  with  KPP  and  TP  (Figs.  11  and  12  and  also  Figs,  lib,  e  and 
12b,  e  versus  Fig.  5). 

Fig.  13  shows  the  distribution  of  shear  Richardson  numbers  along  the  flow  direction  averaged 
over  the  spanwise  direction  for  both  the  KPP  and  TP  experiments  at  selected  horizontal  resolution 
of  Ax  =  100  m  and  at  instant  t  -  13,050  s.  The  comparison  is  performed  for  layers  3  and  4  (out  of 
5)  because  Ri  is  not  defined  at  the  top  and  the  bottom  layers  in  TP.  Furthermore,  the  entrainment 
into  layer  2  is  weak.  In  KPP,  Ri  is  defined  at  the  layer  interfaces,  whereas  it  is  defined  at  the  middle 
of  the  layers  in  TP.  In  order  to  remove  this  difference,  a  linear  interpolation  was  used  to  estimate 
Ri  in  the  middle  of  the  layers  in  KPP.  As  shown  in  Fig.  13,  Ri  is  mostly  in  the  range  of  0. 1-0.2, 
which  is  quite  small  compared  to  Ric.  This  finding  explains  why  the  overall  results  are  insensitive 
to  the  value  of  Ric. 

5.4.  Further  examinations 

In  order  to  explore  the  sensitivity  of  our  results  to  mixing  induced  by  the  bottom  shear  stress, 
we  explored  the  impact  of  the  BBL  formulation  given  in  Eq.  (3)  by  running  all  KPP  experiments 
(five  different  resolutions  and  two  different  Kmax)  with  and  without  BBL.  No  difference  has  been 
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E(t)  in  modified  HYCOM-KPP  vs  Nek5000 


E(t)  in  modified  HYCOM-TP  vs  Nek5000 


Fig.  10.  E(t)  in  HYCOM  experiments  with  (a)  modified  KPP  (Km ax  =  2500  cm2  s_1)  and  (b)  modified  TP  (Ca  =  0.15)  at 
different  horizontal  grid  spacings;  (x)  Ax  =  1000  m,  (A)  Ax  =  500  m,  (O)  Ax  =  100  m,  (0)  Ax  =  50  m,  (□)  Ax  =  20  m. 
Entrainment  parameters  from  2D  and  3D  nonhydrostatic  experiments  with  Nek5000  are  shown  in  the  background, 
dotted:  2D,  dashed:  3D. 


found.  Mixing  takes  place  from  above  the  gravity  current,  in  agreement  with  the  picture  put  forth 
by  Peters  et  al.  (in  press)  based  on  observations  of  the  Red  Sea  overflow  in  a  narrow  channel,  in 
which  the  bottom  properties  are  largely  preserved,  and  mixing  is  mostly  confined  to  the  shear 
layer  above  a  thick  and  homogeneous  bottom  layer.  The  BBL  formulation  thus  does  not  play 
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KPP: 


TP: 


(a)  salinity  in  2-D  (X-Z),  t=4950  sec 

0, - 


(b)  t=9450  sec 

0 - 


(c)  t= 13050  sec 

0 - 


x,  km 


(d)  salinity  in  2-D  (X-Z),  t=4950  sec 


Fig.  11.  Evolution  of  salinity  distribution  in  time  in  HYCOM  experiments  with  modified  parameterizations,  (a) 
t  =  4960  s,  (b)  t  =  9450  s,  and  (c)  t  =  13,050  s  using  modified  KPP  (Xmax  =  2500  cm2  s-1)  with  Av  =  1000  m,  and  (d),  (e), 
(f)  using  modified  TP  (CA  =  0.15). 


a  role  in  the  present  set  of  numerical  experiments,  but  it  might  become  more  important  in  over¬ 
flows  subject  to  lateral  spreading. 

We  have  noted  above  that  the  laboratory  experiments  of  Ellison  and  Turner  (1959)  had  rather 
large  slope  angles,  much  larger  than  in  nature.  This  raised  the  question  how  well  the  TP  scheme 
handles  different  slope  angles.  In  order  to  examine  this  issue,  several  additional  experiments  have 
been  conducted  with  slope  angles  of  9  =  2°  and  0=1°  using  modified  KPP  and  TP  formulations  at 
a  selected  horizontal  resolution  of  Av  =  100  m.  The  results  indicate  that  entrainment  curves  E{t) 
from  KPP  remain  virtually  unchanged  in  response  to  changes  in  the  slope  angle  (Fig.  14a),  while 
those  from  TP  exhibit  stronger  variations.  In  response  to  a  3.5-fold  decrease  of  the  slope  angle,  the 
equilibrium  entrainment  parameter  decreases  by  about  20%  in  TP  (Fig.  14b). 


6.  Summary  and  discussion 

Our  understanding  of  the  dynamics  of  overflows  is  based  on  the  results  of  dedicated  observa¬ 
tional  programs  in  the  Mediterranean  Sea  overflow  (Baringer  and  Price,  1997a,b),  Denmark  Strait 
overflow  (Girton  et  al.,  2001,  2003),  Red  Sea  overflow  (Peters  et  al.,  in  press;  Peters  and  Johns,  in 
press),  Faroe  Bank  Channel  (Price,  2004)  and  Antarctic  Ocean  (Gordon  et  al.,  2004),  and  also  of 
laboratory  tank  experiments  (e.g.,  Ellison  and  Turner,  1959;  Simpson,  1987;  Hallworth  et  al., 
1996;  Monaghan  et  al.,  1999;  Baines,  2001;  Cenedese  et  al.,  2004),  and  process  modeling  studies 
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KPP: 


TP: 


(c)  t= 13050  sec 

0 - 


10  15  20  25 

x,  km 


(f)  t=  13050  sec 


10  15  20  25 

x,  km 


Fig.  12.  Evolution  of  salinity  distribution  in  time  in  HYCOM  experiments  with  modified  parameterizations,  (a) 
spanwise-average  at  t-  4960  s,  (b)  spanwise-average  at  t  =  9450  s,  and  (c)  spanwise-average  at  t-  13,050  s  using 
modified  KPP  (Xmax  =  2500  cm2  s-1)  with  Ax  =  20  m,  and  (d),  (e),  (f)  using  modified  TP  (CA  =  0.15). 


(e.g.,  [49-52]).  It  is  important  that  this  knowledge  is  incorporated  in  OGCMs  in  the  form  of 
appropriate  mixing  parameterizations. 

In  this  study,  experiments  are  conducted  using  an  OGCM,  HYCOM,  in  an  idealized  setting 
that  represents  the  basics  of  shear-induced  mixing  in  bottom  gravity  currents,  that  is,  the  flow 
of  a  dense  water  mass  released  at  the  top  of  a  wedge,  which  is  20  km  long,  2  km  wide,  and  has 
a  slope  of  9  =  3.5°  with  respect  to  the  horizontal.  Similar  experiments  have  been  carried  out  by 
dzgokmen  et  al.  (2004a)  using  a  high-order  nonhydrostatic  spectral  element  model,  Nek5000,  a 
general  Navier-Stokes  solver  developed  by  Fischer  (1997).  Our  HYCOM  experiments  are  config¬ 
ured  as  similar  as  possible  to  the  Nek500  setting,  and  are  conducted  with  five  different  horizontal 
resolutions  of  1000  m,  500  m,  100  m,  50  m,  and  20  m.  In  HYCOM,  two  mixing  parameterizations 
are  used,  (i)  KPP  (Large  et  al.,  1994,  1999),  a  class  of  multi-purpose  mixing  algorithms  which  in¬ 
cludes  a  shear-induced  mixing  scheme  based  on  results  from  Pacanowski  and  Philander  (1981), 
and  (ii)  TP,  which  has  been  developed  for  overflows  by  Hallberg  (2000)  based  on  laboratory  re¬ 
sults  from  Ellison  and  Turner  (1959)  and  Turner  (1986).  Both  schemes  are  based  on  the  local  gra¬ 
dient  Richardson  number,  but  they  differ  in  that  a  vertical  diffusivity  is  used  in  KPP  while  the 
entrainment  velocity  is  specified  in  TP.  We  explore  how  results  from  HYCOM  with  KPP  and 
TP  compare  to  each  other  and  to  those  from  Nek5000,  and  whether  results  change  significantly 
as  a  function  of  the  model  resolution. 

It  is  found  that  KPP  results  in  significantly  less  gravity  current  entrainment  than  that  in  the  ref¬ 
erence  experiment  with  Nek5000,  while  TP  leads  to  significantly  more  entrainment  than  both. 
Specifically,  the  entrainment  parameter  (defined  in  Section  5.2)  converges  to  £«lx  10-3  in 
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KPP,  Layer  #3  TP,  Layer  #3 


KPP,  Layer  #4  TP,  Layer  #4 


Fig.  13.  The  distribution  of  Ri  along  the  flow  direction  averaged  over  the  spanwise  direction  in  the  case  of  KPP  in 
layers  (a)  3  and  (c)  4,  and  in  the  case  of  TP  in  layers  (b)  3  and  (d)  4,  at  the  selected  horizontal  resolution  of  Av  =  100  m 
and  at  t  -  13,050  s. 

experiments  with  KPP,  E  «  8  x  1 0  in  experiments  with  TP,  whereas  E  ~4x  10-3  in  the  3D  exper¬ 
iment  with  Nek5000.  The  results  are  fairly  independent  of  the  horizontal  grid  resolution.  KPP  and 
TP  are  then  tuned  using  results  from  Nek5000,  and  it  is  found  that  this  requires  an  increase  of  Kmax 
from  50  cm2  s'1  to  2500  cm2  s~  1  for  KPP,  and  a  decrease  of  Ca  from  1  to  0.15  for  TP. 

Given  that  the  parameters  of  KPP  and  TP  needed  to  be  changed  significantly  in  order  for  them 
to  match  the  results  from  high-resolution  nonhydrostatic  3D  model  runs,  further  discussion  of  the 
structure  of  these  parameterization  schemes  is  needed.  With  respect  to  the  original  experiments 
analyzed  by  Ellison  and  Turner  (1959)  and  Turner  (1986),  which  underlie  TP,  one  could  raise 
questions  concerning  viscous  effects  of  bottom  and  side  boundaries  in  a  very  small  tank  2  m  long 
and  10  cm  wide,  and  a  flow  only  10  cm  deep.  Other  questions  concern  the  large  slopes  in  the  tank 
experiments,  12°  ^  9  ^  90°,  compared  to  small  angles,  0°  <  5°  in  nature.  However,  the  fundamen¬ 
tal  reason  why  the  parameter  values  of  Hallberg’s  (2000)  TP  scheme  have  to  be  considered  adjust¬ 
able  is  that  they  were  taken  unchanged  from  an  algorithm  employing  the  bulk  Richardson 
number  of  a  single-layer  bottom  gravity  flow  and  applied  in  a  new  algorithm  employing  gradient 
Richardson  numbers  in  a  multi-layered  shear  flow.  The  differences  between  these  two  physical  set¬ 
tings  are  vast. 

Numerical  experiments  by  Papadakis  et  al.  (2003)  provide  further  incentive  to  tune  TP.  They 
conducted  simulations  of  the  Mediterranean  Sea  outflow  using  HYCOM  with  TP  with  encourag¬ 
ing  results.  However,  in  order  to  obtain  a  realistic  path  of  the  overflow  and  to  achieve  the  gene¬ 
ration  of  subsurface  eddies  (Meddies),  in  a  somewhat  ad  hoc  approach  they  resorted  to  applying 
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E(t)  in  modified  HYCOM-KPP  for  0  =  1°,2°,3.5° 


E(t)  in  modified  HYCOM-TP  for  0  =  10,20,3.5° 


Fig.  14.  E{t)  in  HYCOM  experiments  with  (a)  modified  KPP  (Kmax  =  2500  cm2  s_1)  and  (b)  modified  TP  (CA  =  0.15)  at 
Ax  =  100  m  for  different  slope  angles  9=  1°  (black  lines),  9-  2°  (blue  lines),  9  =  3.5°  (red  lines).  Entrainment 
parameters  from  2D  and  3D  nonhydro  static  experiments  with  Nek5000  are  shown  in  the  background,  dotted:  2D, 
dashed:  3D.  (For  interpretation  of  the  references  in  colour  in  this  figure  legend,  the  reader  is  referred  to  the  web  version 
of  this  article.) 

the  TP  scheme  only  every  144th  time  step  rather  than  at  every  step.  Further,  when  TP  is  used  as  a 
general  shear-induced  mixing  parameterization  in  North  Atlantic  simulations  of  HYCOM,  it 
leads  to  unrealistically  high  mixing  rates  in  the  equatorial  regions. 
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Problems  with  KPP  have  also  become  obvious.  The  KPP-modeled  Mediterranean  out-flow 
sinks  far  deeper  than  observations  in  recent  high-resolution  (1/12°)  simulations  using  HYCOM, 
indicating  that  mixing  induced  by  KPP  is  insufficient  when  directly  applied  to  overflows.  We  have 
already  noted  that  KPP  cannot  be  universally  valid  because  of  its  simplistic  structure,  a  dimen¬ 
sional  constant  times  a  function  of  a  nondimensional  parameter,  Ri.  Specifically,  Kmax  was  deter¬ 
mined  from  LES  modeling  of  the  diurnal  cycle  of  surface  mixed  layer  at  the  equator  subject  to  a 
specific  forcing.  Physical  intutition  leads  to  the  expectation  that  Kmax  should  vary  with  the 
strength  of  the  forcing  and  that  it  should  not  be  expected  that  this  particular  value  of  Kmax  to  hold 
in  bottom  gravity  current  mixing. 

Noting  that  the  original  development  of  KPP  (Large  and  Gent,  1999)  in  addition  to  LES  sim¬ 
ulations  also  contemplated  the  oceanic  turbulence  observations  of  Peters  and  Gregg  (1988),  we 
reviewed  their  measurements  in  the  light  of  our  current  study.  Within  the  high-shear,  low-7?/  set¬ 
ting  of  the  Pacific  Equatorial  Undercurrent  at  140°  W,  the  1984  Tropic  Heat  1  Experiment  found 
much  stronger  mixing  at  night,  when  the  ocean  loses  heat  and  the  surface  mixed  layer  undergoes 
convection,  than  during  daytime,  when  the  solar  heat  input  stabilizes  the  upper  ocean.  Hence,  in 
this  environment  the  forcing  of  the  turbulence  has  a  nighttime  maximum.  Fig.  15  depicts  hourly 
averages  of  the  eddy  diffusivity  of  heat,  Kh,  as  function  of  the  local  Ri  separately  for  daytime  and 
nighttime.  While  the  overall  shape  of  the  average  curve  Kh  =  Kh(Ri )  does  not  change  significantly 
between  day  and  night,  nighttime  adds  large  Kh  to  the  high-Kh  end  of  the  curve.  This  is  like  vary¬ 
ing  Kmax  in  KPP.  Therefore,  turbulence  parameterizations  should  include  both  a  dependence  on 
the  forcing  and  a  dependence  on  the  flow  Richardson  number.  This  requirement  holds  for  TP  but 
not  for  KPP. 
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Fig.  15.  Hourly  averages  from  the  1984  Tropic  Heat  1  Experiment  in  the  Pacific  Equatorial  Undercurrent  at  0°,  140° 
W,  eddy  diffusivity  of  heat  versus  arc  tangent  of  the  local  gradient  Richardson  number,  (a)  daytime  data  subject  to 
oceanic  heat  gain  and  upper  ocean  stabilization,  (b)  nighttime  data  subject  to  oceanic  heat  loss  and  mixed  layer 
convection.  The  data  span  the  upper  ~150  m  depth  with  the  core  of  the  Undercurrent  and  a  minimum  in  shear  near 
110  m.  Note  that  the  shape  of  the  data  scatter  changes  little  between  daytime  and  nighttime,  while  more  large  Kh  appear 
at  the  high-Zh  end  of  the  curve  at  night.  The  vertical  lines  in  (a)  and  (b)  indicate  Ri  =  1/4. 
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The  preceding  requirement  is  consistent  with  common  and  accepted  turbulence  parameteriza- 
tions  more  complex  than  KPP  and  TP.  Two-equation  turbulence  closures  of  all  varieties  (e.g., 
Mellor  and  Yamada,  1982;  Baumert  and  Peters,  2004)  represent  the  /^/'-dependence  as  “stability 
functions,”  while  the  dependence  on  the  forcing  is  handled  by  the  pair  of  predictive  differential 
equations  for  the  turbulent  kinetic  energy  and  another  variable  related  to  the  turbulent  length 
scale.  Sub-grid  scale  schemes  commonly  employed  in  LES  models  are  similar,  if  simpler.  Follow¬ 
ing,  e.g.,  Smagorinksy  (1963),  Deardorflf  (1973),  Schumann  (1991)  and  Stevens  et  al.  (2000),  one 
can  write 


K  =  cl2  \  S\f(Ri),  (14) 

where  S 2  =  f^A/  (i,j  =  1,2,3)  is  the  resolved  strain  rate.  Further,  Ay  -  dujdxj  +  dujldxi  denotes 
the  resolved  scale  deformation  using  a  grid  spacing  proportional  to  /  (typically,  /  =  (AvAyAz)173), 
and  c  is  an  empirical  constant.  The  effect  of  stratification  is  incorporated  by  specifying  a  monot¬ 
onously  decreasing  function  that  satisfies 


rm 


1  for  Ri  =  0, 

0  for  Ri  ^  Ric. 


(15) 


The  detailed  shape  of  this  curve  would  require  additional  information  about  the  mixing  process 
(e.g.  as  in  the  Peters  and  Gregg  (1988)  observations),  but  even  a  linear  relationship  could  suffice  as 
a  first-order  approximation. 

The  key  point  is  that  the  dynamical  factor  determining  the  amplitude  of  the  mixing  coefficient 
at  low  Ri,  which  is  a  function  of  the  resolved  strain  rate,  l2\S\,  is  replaced  by  a  peak  diffusivity  of 
^|i?/=o  -  Anax  -  50  cm2  s-1  in  KPP,  based  on  results  from  a  physical  regime  quite  different  than 
oceanic  overflows.  In  contrast,  by  relating  wE  to  A  U,  TP  does  avoid  a  hard  limit  for  peak  effective 
diffusivity,  and  the  implied  diffusivity  includes  both  a  dependence  on  the  flow  Richardson  number 
and  on  the  forcing  via  the  resolved  model  local  velocity  structure  A  U.  This  explains  why  the 
extrainment  parameter  in  TP  changes  in  response  to  variations  in  the  slope  angle  whereas  KPP 
does  not  seem  to  show  any  response. 

In  future  studies,  it  will  be  explored  how  KPP,  a  popular  mixing  model  for  many  OGCMs,  can 
be  modified  to  incorporate  the  dependence  of  mixing  coefficients  on  the  forcing  for  overflow  sce¬ 
narios.  Also,  it  needs  to  be  further  investigated  how  accurately  TP  captures  the  dependence  of  the 
entrainment  on  the  slope  angle. 

Finally,  the  environment  used  in  this  study — as  well  as  in  the  experiments  of  Ellison  and  Turner 
(1959) — is  homogeneous,  whereas  ambient  stratification  can  have  a  significant  effect  on  the  nature 
of  the  mixing  process  and  entrainment  (e.g.,  Baines,  2001),  and  therefore  the  dynamics  of  over¬ 
flows.  In  order  for  mixing  parameterizations  to  be  applicable  to  the  ocean,  the  effect  of  ambient 
stratification  needs  to  be  taken  into  account.  This  issue  will  also  be  explored  in  the  near  future. 
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